# Set working directory to whatever directory the data is in
# Example:
# setwd("~/Dropbox/alliance_reputation/replication")

# make sure you have all of these packages installed.
# use command 
# install.packages(c("ggplot2", "dplyr", "arm", "foreign", "games", "sandwich", "lmtest", "MASS", "Zelig", "grid"), dependencies=T)

# clear memory in R
rm(list=ls())

# load libraries
library(ggplot2)
library(dplyr)
library(arm)
library(foreign)
library(games)
library(sandwich)
library(lmtest)
library(MASS)

#read in data
alliance <- read.dta("alliance_replication.dta")

############################################
# Create Figure 1 from the main manuscript #
############################################

fig1 <- ggplot(data= subset(alliance, regyrvio_avg > 0), aes(x=year, y=regyrvio_avg))
fig1 <- fig1 + geom_point() + geom_smooth()
fig1 <- fig1 + scale_x_continuous(breaks=seq(1920, 2000, by=20))
fig1 <- fig1 + facet_grid( . ~ region1)
fig1 <- fig1 + theme_bw() + xlab("Year") + ylab("Percent of Countries Ending \n an Alliance Through Violation") 

ggsave(fig1, file="Figure_1.pdf", width=10, height = 3)


############################################
# Create Figure 2 from the main manuscript #
############################################

# function to calculate clustered standard errors
# same as the "cluster"" command in stata
clx <- function(fm, dfcw, cluster){
  M <- length(unique(cluster))
  N <- length(cluster)
  dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
  u <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
  vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
  coeftest(fm, vcovCL) 
}

# clx2 returns a variance covariance matrix for cluster robust se's
clx2 <- function(fm, dfcw, cluster){
  M <- length(unique(cluster))
  N <- length(cluster)
  dfc <- (M/(M-1))*((N-1)/(N-fm$rank))
  u <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum))
  vcovCL <- dfc*sandwich(fm, meat=crossprod(u)/N)*dfcw
  vcovCL
}


# Fit the 6 model's from Figure 2 in the main manuscript (from left to right)
# The coeffiecients and standard errors from these models will be used to create Figure 2.
# The standardize command rescales the regression inputs so that they 
# are centered at 0 and normalized by 2 sd's
# The coeffients for these normalized inputs form the basis for our 
# graph "Figure 2" in the main manuscript.

# Create Matrix to store coefficients
coef.matrix <- matrix(nrow = 6, ncol = 6)

# Create Matrix to store clustered standar errors
se.matrix <- matrix(nrow = 6, ncol = 6)

# Model 1 and 2 comparing our measures with state clustered SE
# but no time dynamics

#Model 1
syrfe.cumvio.notime <- standardize(glm(alliance_gain ~  cum_vio_dum + majpow + democ_dum + lag_num + I(lag_num^2) + as.factor(region1), family= binomial(link=logit),  data=alliance, x=TRUE))

#estimated coefficients and se's using cluster robust se's assuming correlation
# within state observations
clustse.syrfe.cumvio.notime <- clx(syrfe.cumvio.notime, 1, alliance$state[!is.na(alliance$lag_num)])

# store model coeffients and clustered se's 
coef.matrix[2:6,1] <- clustse.syrfe.cumvio.notime[2:6,1]
se.matrix[2:6,1] <- clustse.syrfe.cumvio.notime[2:6,2]


#variance covariance matrix
vcov.clustse.syrfe.cumvio.notime <- clx2(syrfe.cumvio.notime, 1, alliance$state[!is.na(alliance$lag_num)])


#Model 2
syrfe.rep.notime <- standardize(glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num + I(lag_num^2) + as.factor(region1), family= binomial(link=logit),  data=alliance, x=TRUE))

clustse.syrfe.rep.notime <- clx(syrfe.rep.notime, 1, alliance$state[!is.na(alliance$lag_num)])

# store model coeffients and clustered se's 
coef.matrix[c(1,3:6),4] <- clustse.syrfe.rep.notime[2:6,1]
se.matrix[c(1,3:6),4] <- clustse.syrfe.rep.notime[2:6,2]

### Perform Clarke test for 
### non-nested models on Model's 1 and 2 as described on
### page 29 of the manuscript
clarke(syrfe.cumvio.notime, syrfe.rep.notime)


# Model 3 and 4 comparing our measures with state clustered SE
# plus time dynamics

#Model 3
syrfe.cumvio.time <- standardize(glm(alliance_gain ~  cum_vio_dum+ majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + as.factor(region1), family= binomial(link=logit),  data=alliance))

clustse.syrfe.cumvio.time <- clx(syrfe.cumvio.time, 1, alliance$state[!is.na(alliance$lag_num)])

# store model coeffients and clustered se's 
coef.matrix[2:6,2] <- clustse.syrfe.cumvio.time[2:6,1]
se.matrix[2:6,2] <- clustse.syrfe.cumvio.time[2:6,2]



#Model 4
syrfe.rep.time <- standardize(glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + as.factor(region1), family= binomial(link=logit),  data=alliance))

clustse.syrfe.rep.time <- clx(syrfe.rep.time, 1, alliance$state[!is.na(alliance$lag_num)])

# store model coeffients and clustered se's 
coef.matrix[c(1,3:6),5] <- clustse.syrfe.rep.time[2:6,1]
se.matrix[c(1,3:6),5] <- clustse.syrfe.rep.time[2:6,2]


#Model 5 and 6  add a decay term

# Model 5
syrfe.cumvio.decay <- standardize(glm(alliance_gain ~  cum_vio_dum + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + decay_gen + I(decay_gen^2) + I(decay_gen^3)+ as.factor(region1), family= binomial(link=logit),  data=alliance))

clustse.syrfe.cumvio.decay <- clx(syrfe.cumvio.decay, 1, alliance$state[!is.na(alliance$lag_num)])

# store model coeffients and clustered se's 
coef.matrix[2:6,3] <- clustse.syrfe.cumvio.decay[2:6,1]
se.matrix[2:6,3] <- clustse.syrfe.cumvio.decay[2:6,2]


# Model 6
# unnormalized model for later
syrfe.rep.decay <- glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + decay_gen + I(decay_gen^2) + I(decay_gen^3)+ as.factor(region1), family= binomial(link=logit),  data=alliance)
clustse.syrfe.rep.decay <- clx(syrfe.rep.decay, 1, alliance$state[!is.na(alliance$lag_num)])

# normalized for Figure2 
syrfe.rep.decay.norm <- standardize(glm(alliance_gain ~  rep_dum_mem_indef + majpow + democ_dum + lag_num  + I(lag_num^2)  + t + I(t^2) + I(t^3)  + decay_gen + I(decay_gen^2) + I(decay_gen^3)+ as.factor(region1), family= binomial(link=logit),  data=alliance))
clustse.syrfe.rep.decay.norm <- clx(syrfe.rep.decay.norm, 1, alliance$state[!is.na(alliance$lag_num)])

# store model coeffients and clustered se's 
coef.matrix[c(1,3:6),6] <- clustse.syrfe.rep.decay.norm[2:6,1]
se.matrix[c(1,3:6),6] <- clustse.syrfe.rep.decay.norm[2:6,2]

# Run additional code to make  Figure 2 
source("table_graph_code.R")


######################################
# Create Figure 3 from the main text #
###################################### 
# This graphs the estimated effect of reputation
attach(alliance)
sims <- 5000
Bhat <- as.vector(clustse.syrfe.rep.decay[,1])
par.draws <- mvrnorm(sims, mu=Bhat, Sigma=clx2(syrfe.rep.decay, 1, alliance$state[!is.na(alliance$lag_num)]))

# set all variables to their median
X1 <- cbind(1, 0, median(majpow, na.rm=T), median(democ_dum, na.rm=T), median(lag_num, na.rm=T), median(lag_num, na.rm=T)^2, median(t, na.rm=T), median(t, na.rm=T)^2, median(t, na.rm=T)^3, median(decay_gen, na.rm=T), median(decay_gen, na.rm=T)^2, median(decay_gen, na.rm=T)^3, 0,0,0,0)

est.fx <- c()
for(i in 6:81){
  X1[2] <- i
  mu <- 1/(1+exp(-(par.draws%*%t(X1))))
  mean.fx <- mean(mu)
  sd.fx <- sd(mu)
  ci.fx <- quantile(mu, c(.025, 0.975))
  est.fx <- rbind(est.fx, c(mean.fx,sd.fx,ci.fx))
}
est.fx <- as.data.frame(est.fx)
names(est.fx) <- c("mean","sd","ci.lo", "ci.hi")

fx.plot <- ggplot(data=est.fx, aes(x=seq(6,81,1), y=mean))+
	geom_ribbon(aes(x=seq(6,81,1), ymin=ci.lo, ymax = ci.hi, alpha=.3))+
	geom_line(size = 2)+
	xlab("Unreliability Score")+ylab("Pr(Alliance Gain | 1 Violation)")

fx.plot <- fx.plot + ggtitle("The Probability of Gaining a new Alliance\nas a Function of Unreliability") +theme_bw() +  theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank(), legend.position="none")

ggsave(fx.plot, file="Figure_3.pdf", width= 6, height =6)

####################################
# Create figure A1 in the appendix #
####################################

timeseq <- 0:max(alliance$t, na.rm=T)

# Set all variables to their median
newX2 <- cbind(1, median(rep_dum_mem_indef, na.rm=T), median(majpow, na.rm=T), median(democ_dum, na.rm=T), median(lag_num, na.rm=T), median(lag_num, na.rm=T)^2, timeseq, timeseq^2,timeseq^3,0,0,0, 0,0,0,0)
Bhat2 <- as.vector(clustse.syrfe.rep.decay[,1])
XBhat2 <- as.matrix(newX2)%*%Bhat2
prob2 <- 1/(1+exp(-XBhat2))

time.plot <- qplot(x=timeseq, y=prob2, geom="line", xlab="t", ylab="Pr(Alliance Gain | t, median(X)")

time.plot <- time.plot + ggtitle("Probability of Gaining a New Alliance\nAfter an Alliance Gain") +theme_bw() +theme( panel.grid.minor = element_blank(), panel.grid.major = element_blank())

ggsave(time.plot, file="Figure_A1.pdf", height=6, width=6)

####################################
# Create Figure A2 in the appendix #
####################################

#graph decay
decayseq <- 0:max(alliance$decay_gen, na.rm=T)

#Set all variables to their median
newX <- cbind(1, median(rep_dum_mem_indef, na.rm=T), median(majpow, na.rm=T), median(democ_dum, na.rm=T), median(lag_num, na.rm=T), median(lag_num, na.rm=T)^2,median(t, na.rm=T), median(t, na.rm=T)^2, median(t, na.rm=T)^3, decayseq, decayseq^2,decayseq^3, 0,0,0,0)
Bhat <- as.vector(clustse.syrfe.rep.decay[,1])
XBhat <- as.matrix(newX)%*%Bhat
prob <- 1/(1+exp(-XBhat))

decay.plot <- qplot(x=decayseq, y=prob, geom="line", xlab="t", ylab="Pr(Alliance Gain | t, median(X)")
decay.plot <- decay.plot + ggtitle("Probability of Gaining a New Alliance\nAfter an Alliance Violation")+theme_bw()+theme(panel.grid.minor = element_blank(), panel.grid.major = element_blank())

ggsave(decay.plot, file="Figure_A2.pdf", width= 6, height=6)






